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Heat release by controlled continuous-time 
Markov jump processes 



Abstract We derive the equations governing the protocols minimizing the heat 
released by a continuous-time Markov jump process on a one-dimensional count- 
able state space during a transition between assigned initial and final probability 
distributions in a finite time horizon. In particular, we identify the hypotheses on 
the transition rates under which the optimal control strategy and the probability 
distribution of the Markov jump problem obey a system of differential equations 
of Hamilton-Jacobi-Bellman-type. As the state-space mesh tends to zero, these 
equations converge to those satisfied by the diffusion process minimizing the heat 
released in the Langevin formulation of the same problem. We also show that in 
full analogy with the continuum case, heat minimization is equivalent to entropy 
production minimization. Thus, our results may be interpreted as a refined version 
of the second law of thermodynamics. 



1 Introduction 

Molecular motors and more generally nano-machines operate in viscous, fluctuat- 
ing environments. It is therefore useful, if not necessary, to model these systems by 
means of stochastic processes and to describe their behavior in terms of probabil- 
ity distributions. Fluctuation relations I.15..16..18„24„31-!32::10j22J impose general 

Paolo Muratore-Ginanneschi 

University of Helsinki, Department of Mathematics and Statistics, P.O. Box 68 FIN-00014, 

Helsinki, Finland. E-mail: paolo.muratore-ginanneschi@helsinki.fi 

Carlos Mejia-Monasterio 

Laboratory of Physical Properties, Department of Rural Engineering, Technical University of 

Madrid, Av. Complutense s/n, 28040 Madrid, Spain. E-mail: carlos.mejia@upm.es 

Luca Peliti 

Dipartimento di Scienze Fisiche and Sezione INFN, Universita Federico II, Complesso Monte 

S. Angelo, 180126 Napoli, Italy. E-mail: peliti@na.infn.it 



constraints on these probability distributions which can be and have been exten- 
sively tested experimentally see e.g. Il46ll33l l6ll43ll9l l20l . A unified approach to 
these results from a theoretical standpoint can be found in ff] while a review with 
emphasis on experimental aspects and extensive reference to the literature is [40 1. 
Both in experimental and theoretical investigations of nano-machines it is crucial 
to distinguish between fast fluctuating configurational variables and control pa- 
rameters, i.e., variables whose state is determined by extemal macroscopic sources. 
For example, in an experiment following the trajectory of a micron-sized bead im- 
mersed in water and captured in an optical trap, the control parameter could be the 
center of the trap measured in the laboratory frame whereas the configurational 
variable is the displacement of the bead P31 [T1. 

An important observation made in 0411 is that the knowledge of the control pa- 
rameters driving a nano-system in a finite-time transition between two assigned 
states while minimizing a suitable non-equilibrium thermodynamical functional 
(the mean work done on the system in the examples considered by 14111 '). yields 
substantial improvements in the measurement from finite-time path sampling of 
thermodynamical indicators (e.g. free energy differences) in both numerical and 
experimental studies. 

Recent results fTT'2l have shown that the problem posed in [41] admits a precise 
mathematical formulation in the language of optimal control theory see e.g. 1 17) 
and also |44| for a short, brilliant introduction. A conspicuous consequence is that 
for fast processes described by a Langevin dynamics the minimization of the heat 
released or of the work done in a transition between assigned states can be exactly 
mapped into Monge-Ampere-Kantorovich optimal transport problems 1451 . 
The aim of the present contribution is to extend these results to transitions de- 
scribed by continuous time Markov jump processes on a countable state space 1127 1 
[26 1 . Markov jump processes have been often considered in the study of fluctua- 
tion theorems see e.g. 132, 23 , 34 , 13 , 141 as model problems for non-equilibrium 
thermodynamics. The application closest to the scope of the present contribution 
can be found in 1 12|. Namely, the authors of fT?! modeled a single-level fermion 
system interacting with a thermal reservoir by a two-state Glauber jump-process 
with the aim of determining the protocol raising the energy level with minimal 
work done on the system. Rigorous optimal control theory for jump processes has 
been developed long ago II39II11II . Adapting it to the minimization of the heat re- 
leased during a transition between states leads to a formulation in weak sense of 
the variational problem similarly to what happens in stochastic mechanics 1211 
and Euclidean quantum mechanics II48L Heat release minimization is, within our 
working hypotheses, equivalent to entropy production minimization. Thus, our re- 
sults have a straightforward physical interpretation as refined bounds for the sec- 
ond law of thermodynamics [4 2] amenable to direct experimental testing [5[. In 
this respect, the aforementioned distinction between configurational and control 
parameters has an immediate, important consequence for control theory. Identify- 
ing without any further specification the jump rates of the process as the controls 
implies in physical terms acting at the fastest possible time scale of the system. 
An intuitive and, somewhat trivial, consequence is that optimal control will be 
a jump process. In our view, a more physically relevant approach is to instead 
inquire entropy production minimization with respect to the broadest set of pro- 
cess parameters which lead to smooth, macroscopic optimal control protocols. 



The mathematical a-priori condition for this distinction is coercivity, a well known 
concept in control theory IITtI : the convexity of the cost functional with respect to 
the control. Non-coercivity, e.g. linearity of the cost functional with respect to an 
optimization parameter brings generically about singular controls. We interpret 
here singular controls as optimization protocols acting on configurational param- 
eters of the system and, as such of questionable physical realizability. We will 
argue in the conclusions of this paper that the introduction of a control time scale 
originating from the distinction between configurational and control parameters is 
not a peculiarity of Markov jump process but it is also implied in the Langevin 
dynamics formulation |l3]|4]|2l- 

The structure of the paper is as follows. In § |2] we shortly recall basic proper- 
ties of Markov jump processes and introduce the heat functional as a measure of 
irreversibility in transition between states. In order to simplify the notation, we 
restrict the attention to a one-dimensional state space. As noticed in |4|, the heat 
optimal control problem is most naturally formulated in the current velocity for- 
malism [38] see also fS" 19] for application to fluctuation relations. In §[3]we state 
the optimization problem under the hypothesis that transition rates satisfy a local 
detailed balance. We show that in such a case the heat is a convex functional of 
the control only if a certain statistical indicator we call "reduced traffic" fi% is 
bounded. Under this further hypothesis the optimal control satisfy a system of dif- 
ferential equations of Hamilton-Jacobi-Bellman type |17|. In §|4]we derive these 
equations which we show in §|5]to recover in the limit of vanishing lattice spacing 
the optimal mass transport equations of ||3]|4]|2]. Finally, in §|6]we consider the 
two and three state dynamics. For the two state dynamics we show that our opti- 
mal protocol corresponds to an entropy production lower than the one generated 
by optimizing the Glauber jump process as in III21 . This is not surprising because 
using Glauber transition rates is equivalent to fixing the value of the reduced traf- 
fic. We turn then to the case of a three-state space for which we solve numerically 
the optimal transport equations. 



2 Continuous time Markov jump process 

Let ^t a continuous-time Markov jump process taking values on a one-dimensional 
countable state space §. Let /(x, t) any test function 

/:§xM+i-^M (1) 

differentiable at least once with respect to the time variable t G IR+. The mean 
forward derivative of / along the realizations of ^t specifies the generator of the 
process 



iTo'^^-H dt :=(5./ + L/)(x,i) 



dtfix, i) + ^ [/(X, t) - /(x, t)] K,(i|x) (2) 



ies 



The jump rates K^'s are positive definite ( Kf(x|x) > for any pair of states 
X, x) and vanishing on the diagonal ( Kt(x|x) = for all x e §). In particular. 



Kj (x|x) denotes the jump rate from the state x to x. We allow here the jump rates to 
depend upon the time t. The knowledge of the generator characterizes the process. 
Namely, the evolution of the probability distribution 

P(6 = x):=m(x,i) (3) 

is govemed by the Master equation 

-^ (x, t)^J2 {Kt(x|x) m (x, t) - Kt(x|x) m (x, t)} (4) 

which for assigned Kj's is a system of non-autonomous differential equations, 
infinite dimensional if S comprises an infinite number of states. We require the 
transition rates to be sufficiently regular for (S) to admit a unique solution for any 
initial probability distribution TOo(x) assigned at time t = to- It is then expedient 
to introduce the transition probability Pt,ta ^s the semi-group solving 



dPt.t^ , 
dt ^ 


x|Xo) : 


= )^{Kt(x|x)Pt,t„(x|x„) 

xes 

limPi,t„ (x|xo) = 1 



(5b) 



and to describe the evolution of probability measures as 



. (x, t)^Yl Pt^t^ (x|i) m„(X) := (Pt,t„m„)(x) (6) 

xGS 



2. 1 Time-reversal of the probability measure 

On any finite time horizon T = tf —to, under rather general regularity assumptions 
on the transition rates Kt's to any Markov jump process ^ = {^t,t G [to,tf]} 
evolving from an initial probability distribution nioix), it is possible to associate a 
time reversed process starting at the time-evolved distribution nif (x) — m{x,tf). 
The corresponding transition probabilities satisfy the time-reversal condition | 



Pt2.ti (X2|xi) TO(xi,ii) = Pt^^t2 (xi|x2) m(x2,t2) (7) 

for any to < ti < t2 < tf. A straightforward calculation using (|7) (sketched in 
appendix |A]l yields the mean backward derivative of the process 



dtiO [ dt J 



:={dt-Lf){x,t) 



= (aj)(x,t)-^[/(x,.)-/(x,t)]M^il:^ (8) 



ies 



whence we identify the transition rates of the time reversed process 



K*(x|x) := 



Kt{x\x)m{x,t) 
m(x, t) 



(9) 



It is propaedeutic to our scopes to consider an ahemative derivation of this classi- 
cal resuh. Let us suppose that the process ^ be adapted to ths forward sub-sigma 
algebra of the natural filtration of a Poisson-clock process 77 = {77*, t £ [to, i/]} 
statistically invariant under time-reversal and specified by a spatially uniform 
jump-rate r^. Under these hypotheses, Girsanov formula (see e.g. r26"291 or for- 
mula (2.5) of 1 34 1) provides us with two equivalent representations of the Radon- 
Nikodym derivative of the measure V^ of ^ with respect to that V,^ of 77. On the 
one hand, we can write dV^/dVjj in the form of an T'^ -martingale with respect to 
the time to 



dV^ 
dVr, 



iv) 



exp ■ 



dt 



Y.M^\vt)~n 



.xes 



^ \nKtivt\m-)\ (10) 

t6JI(r;) 



On the other hand, Girsanov formula yields also the expression for the Radon- 
Nykodim derivative of any process ( = {Q,t ^ [to,tf]} and transition rates Z, 
absolutely continuous with respect to 77 and, adapted to the its backward filtration 
in the form of an T-",, -martingale with respect to the time tf 



dVc 
dVr, 



iv) 



exp 



dt 



^Zt(x|77t) -Tt 



xes 



Y, lnZ,(7;,-|7/,)^ (11) 

teJ(r)) 



In writing jlOb . (Illb we adhere to the convention of considering cddlag (equiva- 
lently, corlol: continuous on (the) right, limit on (the) left) paths. Furthermore, we 
denote by 1(77) the path-dependent, at most countable set of jumps encountered by 
the realizations of tj. The requirement 



^ law J. 



(12) 



translates then into the condition 



mo('7t„) -J^iv) = ^fiVtf) -j^iv) 



dVn 



dV„ 



(13) 



Upon recalling that increments of a smooth function / evaluated along a "corlol" 
step-path, such as the realizations of a pure jump process, are amenable to the 
form 



f{Vt2,t2)- fivt 



,M)= r dtdtf{Tiut)+ Y. if(Vt,t)~f{r,t-,t)] (14) 



teJ())) 



we can solve ( 1131) for the Z's to recover (|9}. Once we determined the analytic ex- 
pression of the time-reversed jump rates, we can use it to construct an auxiliary 
forward process ^ with measure V; absolutely continuous with respect to V^. We 
can then apply again Girsanov formula to evaluate the Kullback-Leibler diver- 
gence II30II of V; from V^ when the two processes start from the same initial data 
at time to'. 

E«) I - r dt J2 [K*(x|6) - K*(xl6)] + E ln^#Hl (15) 

The notation E^^^ emphasizes that the average is with respect to the measure V^. 
The Kullback-Leibler divergence ( 1151) provides us with a natural probabilistic in- 
dicator of the asymmetry between the forward and the backward evolution. In 
particular, we will show in the following section that (115b can be identified as the 
entropy production during a non-equilibrium thermodynamic transition from the 
state mo to the state m/ in the time horizon [to,tf]. With this goal in view, we 
observe that a straightforward calculation (see appendix |B]) allows us to cast ( 1151) 
into the form 

ICiV^m) = S{tf) - S{to) + 13 Qt,.t^ (16) 

The first term on the right-hand side coincides with the variation of Gibbs-Shannon 
entropy 

S'(t) = -^(m In m)(x,t) (17) 

xes 

between the states mo and my across the time horizon [to, i/]. The second term is 

P Qt„t. = /*' dt J2 In ^4^K,(x|x) m (x, t) (18) 

If we now interpret, following 02II35I . Qtf,to ^s the heat released by the process 
in the transition between the state mo and nif, and 13 with the inverse of the tem- 
perature in units of the Boltzmann constant, the identification ( fTSl) establishes a 
"bridge relation" between the theory of stochastic processes and non-equilibrium 
thermodynamics (see [34] and references therein). Two observations are in order 
before closing this section. The first is that the mathematical hypothesis of abso- 
lute continuity guarantees that the ratios y^t{£,t\£,t-) /^t{£,t\£,t-) in ( 115b are well 
defined. Physically we may identify this condition with that of "local detailed 
balance" 0411 . The second observation is that the relation between the Kullback- 
Leibler divergence (I15t and the entropy production is not limited to Markov jump 
processes but admits a straightforward extension to diffusion processes 1,36 J . 



3 Heat release as a "cost functional" 

Let us consider now two physical states described by assigned probability distri- 
butions nio and nif. We are interested in determining the rates Kj's driving the 
transition between nio and nif in a fixed and finite time horizon T — tf — to such 
that the heat released in the process has a minimum. Physical intuition requires 
that the problem be well-posed if ( 118b provides a good definition of the thermo- 
dynamical heat. This means that the heat must be, as a functional of the transition 
rates, bounded from below so to specify a well-defined "cost function" for the 
aforementioned optimal control problem ifTTI . This is indeed the case by virtue of 
a result of (M\. Using probability conservation 

^(Ltm)(x,t)-0 (19) 

it is possible to couch the integrand in dlSI l into the form 



Y: ln^4^K,(x|x)™(x,t) = a{t) ^{t) (20) 

^^ Kt(x|xj at 



where 



ait) , 

E 



x.xGS 



Kt(x|x)m(x,t) - Kt{x\x)m{x,t) ^^ Kt{x\x)m{x,t) ^ ^ 



Kt(x|x)r7i(x,t) 



can be identified as the entropy production rate. Note that (121b can be regarded as 
a consequence of the positive definiteness of the Kullback-Leibler divergence (115b 
with which it coincides. The integral version of (120b 



n. 
or equivalently the identification 



Q = I ' dt' a{t') - [Sitf) - Sito)] (22) 



STot.it f) - STotito) := /3 Q + S{tf) ~ S{to) = / ' dt'a{t') > (23) 

has two important consequences: 

i if we interpret (3 Qtf.t^ as the entropy variation at constant temperature /3^^ of 
the environment during the transformation, (|2J| is the expression of the second 
law of thermodynamics p2H 7ll34l: 

ii for transition between given states, optimal heat control |l4l|2]| reduces effec- 
tively to entropy production minimization as the boundary conditions fully 
specify the variation of the Gibbs-Shannon entropy. 



As the logarithm of a strictly positive definite matrix can always be expressed 
as the sum of a symmetric and an antisymmetric matrix, we write the transition 
rates as 

Kt(x|i) = G(x,x,i)e^^^ (24) 

where for any fixed time t,G: SxSi— ^Misa positive definite, symmetric function 
of X, X, such that G (x, x,t) = Q for any x and the function A is antisymmetric in 

X, X. 

If we represent the probability distribution of the Markov jump process into 
the form 

m(x,t) = e-^(^'*) (25) 

then 

G(x, x,t) == 7(x,x, i)e 2 (26) 

where 7 using the terminology of 1341 is the "traffic" indicator characterizing fluc- 
tuations far from equilibrium of the jump process. At equilibrium 

m,(x) - e^*(^) , (27) 

the traffic 7 = (m(x)K(x|x) + m(x)K(x|x))/2 measures the total number of jumps 
between x and x. In view of ( l26l) . we will interpret G as the state probability 
distribution discounted component of the traffic and refer to it as "reduced traffic". 
At equilibrium, the antisymmetric function A governs detailed balance relation 
and reduces to 

A,(x,X) = i?4x)-ii',(X) (28) 

Out of equilibrium, A is related to the nonequilibrium driving function F as 

F(x, X, t) := A(x, X, t) - [i?(x, t) - R{x, t)] (29) 

Under these definitions, the relation (|24] |. known as local detailed balance condi- 
tion [25 34], fixes the values of the symmetric jump rates in terms of the equilib- 
rium density and the nonequilibrium driving. 

In terms of reduced traffic and driving function the total entropy variation sat- 
isfies 

S := STot.i^f) ~ StoiX^o) = 

I d< ^ G(x,X,i)F(x,i,t) sinh^^^^i^ [m(i,t)m(x,i)]^/^ (30) 
•'*° x.xes 

We will now show that this representation of the total entropy variation provides 
a well-defined cost functional to describe the minimization of the heat release in 
terms of the control fields respectively specified by the reduced traffic and the 
driving function. In appendix|C]we detail an alternative formulation to the control 
problem closer to the approach followed in [3J. 
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4 Optimal control of the total entropy variation 

Inspection of OOb reveals three important facts. First, the total entropy variation is 
a bilinear form in the probability amplitude 



(j> (x, t) = ^/M^^ (31) 

evolving by (|4]l according to the linear law 

5t0(x,i) = (H0)(x,i) (32a) 



(H 0)(x, t) ■- Y^ G(x, i, t) sinh J^lli^0(i, t) (32b) 



The "Hamiltonian" operator H : § x § x K+ i-^- M is for any fixed time t G M+ 
anti-symmetric under states (x, x) permutation 

Ht = -H (33) 

thus enforcing at any time probability conservation: 

Y,4>\^,t) = l Vt (34) 

As a consequence, the optimal control problem can be treated in full analogy with 
the variational techniques described in [48] in the context of Euclidean quantum 
mechanics. The second fact is that ( 1301) is linear in the reduced traffic. Cost func- 
tional linear in the control are known in general to lead to singular solutions i.e. 
not satisfying smooth (partial) differential equations 0171 . The third fact is that the 
total entropy variation is a convex functional of the current potential F. We expect 
therefore that if the reduced traffic is bounded from below or simply constrained 
to a fixed value, the optimal control problem admits a unique solution with F spec- 
ified by a differential equation of Hamilton- Jacobi-Bellman type [17 "441. Using 
(131b . we write the cost functional specified by the heat release between two given 
states as 

/■*/ ,^ /■*/ 

5:==/ di 2_, '?^(^'OU(x,X,i)0(i,t) EE / di0-U-0 (35a) 



Ffx X t) 
U(x, i, t) := G(x, X, t) F(x, X, t) sinh ,, (35b) 

In ( 135b . we regard the reduced traffic G and the driving function F as independent 
controls only restricted by the requirement that the time boundary conditions on 
the probability amplitudes be satisfied. As a consequence the total variation of the 
cost functional decomposes into 

S' = S'^ + S'^ (36) 
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where 

5i= /'dt{0-U^-0 + 0^-U-0 + 0-U-0^} (37) 

for X = (G, F). In dJ/l i and in what follows we define the first variation as 

oi-i:x'(y,y,«)j;3<|^ (38) 

we also assume that X' has the same parity of X under permutation of its state 
space arguments. Bellman principle states that the optimal Markov control corre- 
sponds to the stationary variation of a local functional J, named the value function, 
of the stochastic process. We will now show that for heat release minimization the 
interpretation of Bellman principle in a weak sense similar to [21] yields the op- 
timal control strategy of physical interest. We start by defining the value function 
as the solution of 

a* J (x, i, f) = [H , J] (x, X, t) - U (x, X, t) (39) 

where as usual 

[H , J] (x, X, f) - ^ {H (x, y, t) J (y, 5, t) - J (x, y, t) H (y, i, i)} (40) 

yes 

We require J to satisfy the final boundary conditions 

J'F(x,x,i/) = jax,i,t„) = (41) 

These conditions stem from the fact that admissible controls are only those driving 
the Markov process between assigned probability distributions at the ends of the 
control horizon. For this reason, we require J(-, i/) to be a pure functional of ^j 
whence ( 141b follows. An immediate consequence of ( 1391 ) is 

5 = (</.-J-0)(fo)-(0-J-0)(t/) (42) 

Using ( [39] l and performing a time-integration by parts we can recast (iJTT i into the 
form 

j-i = (0.j;,.0)(t„) 

+ [' dt4>- (aj - [H , J] + U)^ • = (</. • J^ • 4>){to) (43) 

we can now discriminate between three different cases. 
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4.1 Bellman principle in strong sense 

Interpreted in strong sense, Bellman principle implies 

The stationarity conditions take therefore the general form 

= [H^J]-U;, (45) 

which reduce to the system of equations 

= {5y,xJ(y, i, t) - Jy,,J(x, y, t) - <5y,x(5y,xF(x, X, t)} sinh fiZiM (46a) 



= ^^--^^y' ^' '^ ^^■-^^("' y' ^) G(y, y, t) cosh ^^^^ ^^ ^^ 



5y,x < sir 



jf ;: J • - F(x,x,i) F(x,x,i) F(x,x,t)1 ^ 

-'5y,x%,x <! smh ^ — ^ + -^ — '- cosh ^ — '- \ G(y, y, t) (46b) 



Regarding J as a square matrix in the state-space variables, (I46at imposes that 
only diagonal element be non-vanishing. The condition 

J(i, X, t) - J(x, X, t) - F(x, S, t) = (47) 

is, however, consistent with ( 146b b for G > only if 

F(x,x,t) = Vx,Xe§ (48) 

The conclusion is that an optimal control strategy compatible with the boundary 
conditions must be singular. This means that the entropy production attains its 
infimum on a protocol instantaneously switching at any time during the control 
horizon from a state of equilibrium with the initial condition mo to a state in 
equilibrium with the final m/. We won't delve any deeper here on singular control 
but we will retum to it in the conclusions. 

4.2 Bellman principle in weak sense 
The equality 

(</.-J^-0)(to)- /'rft(0-(U^-[H',J])-0)(t) (49) 

yields the weak-sense stationarity conditions 

O = <^(y,t)F(y,y,i)0(y,i) 

x£S 
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= ^ [(/-(y, i)J(y, X, t)q^ix, t) - (^(x, i)J(x, y, t) 0(y, t)] 



xGS 



>(y, t) 0(y, t) <; 2 tanh fiZlM + F(y, y, t) \ (50b) 



{= 



The conditions are satisfied if J is symmetric under permutation of the state vari- 
ables. Furthermore, if we define the auxiliary field B 



0(y,t)i3(y,t)=^J(y,x,t)0(x,t) (51) 

xGS 

the pair (|50] | simplifies to 

= B(y,t)-B(y,t) + F(y,y,t) (52a) 

= B(y, t) - B{y, t) + 2 tanh fiMlll + F(y, y, t) (52b) 

4.3 Unbounded reduced traffic leads to singular control 

As in the strong sense case, if the reduced traffic can be varied without any lower 
bound ( 1521) is satisfied by ( 1481) . The optimal control strategy is again singular. 

4.4 A constraint on reduced traffic leads to an Hamilton-Jacobi-Bellman 
equations 

A more interesting and presumably more physically relevant situation occurs when 
the reduced traffic is bounded from below by a given constraint. In such a case we 
expect the control to consists of a singular part pushing G to its lower bound and 
then of a smoother part, which minimizes the entropy variation versus F for fixed 
G. The limit case is when G is constrained from the very beginning to a given 
value e.g. 

G(x,X,t) = i-^ (53) 

for tc a constant specifying the characteristic time-scale of the process and the 
variation is taken only with respect to the driving function F. In such a case the 
stationarity condition reduces to (I52bt alone. Upon averaging the equation for the 
value function ( |39l l with respect to a probability amplitude and using (I52bb we 
obtain the equation for the auxiliary field B 

aui ^\ 2 Y^^iZiO • u F(y,x,^) F(y,x,t) 

5,5 (X, t) ^ - >^ ^^ smh ^^— tanh ^^— (54) 

yes ^^ ' ' 
Observing that 

2 F(x,x,i) 



a4B(x,t)-B(x,t)] 



tanh 



atF(x,5,i) (55) 
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we get into the equations satisfied by the driving function 

(9(F(x,x,t) = 

§M1 sinh £^ tanh ^^ - ^ sinh ^^ tanh ^^ 

(56) 



yGS 



1 _ 1 tanh^ ^-^^ 



Degrees of freedom counting reveals, however, that not all transition rates of the 
optimal process can be non- vanishing. Let us suppose first |§| = A^ and then de- 
duce the result for infinite lattice from the limit N ^ oo. Because of probability 
conservation ( 134b . the boundary conditions impose 2N — 2 independent condi- 
tions. By virtue of H^ = — H, the evolution of probability amplitudes is probabil- 
ity preserving and brings forth iV — 1 independent equations. We conclude that 
( 156b can only describe the dynamics of other A^ — 1 degrees of freedom. We iden- 
tify these degrees of freedom by reasoning that A^ — 1 independent non-vanishing 
transition rates are exactly those needed to describe a process jumping only to 
nearest neighbors states. If we posit that the state space lattice spacing is dx > 0, 
we achieve a convenient parametrization of the nearest-neighbor dynamics if we 
define the discrete current velocity 

y(x, i) := F (x + dx, x, t) = -F (x, x + dx, t) (57) 

Gleaning the considerations expounded above we infer that the transition between 
two assigned states releasing the minimal amount of heat is specified by the fol- 
lowing system of Hamilton- Jacobi-Bellman 1 17 1 and probability amplitude trans- 
port equations 



t,dtV{^,t) = 
2 sinh ^ tanh %^ 
2 _ tanh^ ^ 



(j){-x.-\- dx,t) 4>{'X;t) 



4>{yi,t) (j){yL + dx,t) 



2 sinh 


^('^-/-■*)tanh 


2 


(?!>(x — dx,t) 


2 sinh 


2 tanh^ ""^l 
^(^+/^^*) tanh 


1) 

V(yi+dx,t) 
2 


</.(x,i) 
0(x + 2dx,i) 



2_tanh2X(M) </)(x + dx,i) ^^^^^ 



tc5t0(x,i) 



- jsinh Y^^ 0(x + dx, t) - sinh ^^"^ /^""'^^ 0(x - dx, t) \ (58b) 

with boundary conditions 

0(x,t„)=0„(x) & 0(x,i/) = (/./(x) (59) 

The value of the total entropy variation along the stationary control ( I52bb 

S^^2 f — Y^ V{x, t) sinh ^^^' ^' 0(x + dx, t) <j){x, t) (60) 

with the convention (/)(y, i) = if y ^ § in case |§| < cx). 
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4.5 Second variation in the presence of constrained reduced traffic 

We now show that the second variation of the total entropy with respect to the 



driving function is positive definite around the stationary point specified by ( I52bb 
for fixed reduced traffic. For any admissible control, we can write 

J^' = (0 • J'f' • 0) (to) 

= [' dtcp- (U'f' - 2 [H'f , J'f] - [H'f' , J] ) • (61) 

Jta 

The first and third integrands in (I61t give positive definite contributions to the 
second variation when evaluated at stationarity. Namely 

a2u(x,x) 

> (62) 



i9F(x, x,i)9F(x,x,i) 
1 



coshfi^llM + f(flMsinh^("'^'') 



2 4 2 

and, since Hp is anti-symmetric, 

B{x,t) - B{i,t) 



0-[H'f',J]-0- Y, H'F'(x,x,i)0(x,i)0(i,i)- 



x,xGS 



2 



- ^ H'F'(x,x,t)0(x,t)0(X,t) |2tanh^^^^^^ + F(x,x,i)| (63) 

whence it follows 

0-[Hf,J]-0Ix"=o < (64) 

Finally, we observe that 



/ 



dt(0-[H'F,J'F]-0)(i)=O (65) 



since we can re-write it as the trace of the product of a symmetric and an anti- 
symmetric matrix. As a consequence we proved that 

J'f, = f ' dtcl)- (U'f' - [H'f' , J] ) ■ 0|f, > (66) 

■Jt„ 



5 Continuum limit 

The limit of vanishing lattice spacing dx for short distance interactions is most 
straightforward to derive if we posit 

Kt (x ± dx\x) = -^ {2 ± /3 6 (x, i) dx} := M^ (67) 
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We define the continuum limit of the Markov jump process with jump rates iSli 
by letting dx tend to zero while fine-tuning the characteristic time as 



/3 (dxf T 



(68) 



for some finite r > and measuring f3 in units such that /Sdx"^ is non-dimensional. 
In such a case the limit 

(£/)(x):=lim(L/)(x) = 

dxlO 

lim [^ ^^ + ^^) " / (^)] fc+ (^' ^) + [/ (^ - ^^) - / (^)] fe^ (^^ ^) (69) 

dxlo (dx)"^ T 



is finite and yields 






(70) 



which is the generator of the stochastic process described by the stochastic differ- 
ential equation 



,/. X dt /2 , 
d£.t = h{S.t:t) — + J--duJt 



(71) 



driven by a Wiener process Uf 



5.1 Heat density 

Under the hypotheses formulated above, the heat functional reduces to 

Mx\x) 






Kt{x\x) 



Pb{x,t) I3b{x + dx,t) 



2 + b{x,t)dx 
2/3(dx)2r 



(dxY 



(3h{x,t) (3h{x-dx,t) 



2 - b[x, t) dx 
2(3{dx)^T 



0{dxf (72) 



whence finally we get into 



V In 



x^xG 



Kf(a:|x) 



Kt(x|x) 



Pb\x,t)~d^b{x,t) 



+ 0{dx) 



(73) 



We thus recover the expression of the heat density released by Langevin dynamics 
see e.g. fT|. Thus, it is a-priori justified to expect that (158b admits as continuum 
limit the optimal transport equations found in [3^4] • In order to derive explicitly 
this result, we define the continuous limit current velocity v by Taylor expanding 
its lattice counter-part in powers of the mesh dx: 

V(X,t) fix + dx, X.t) , , , 1 „ , X / , X9 ^ . , n3 ,r. .-. 

^ ' ' = — '^-!- = vix, t)dx + -d^v(x, t) idxf + O (dxY (74) 

T T 2 
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5.2 Probability density equation 

Upon inserting ( 174b into (I58bb and rescaling time according to ( 1681) , we obtain 



5t0(x,i) = 

J dxv{x,t) + 2 (^x'v{x,t) 

~ 1 2 (dx)2 T 



{(/)(x,i) +dx(9:j(^(x,i)} 



which yields 

9t(?!)(a;, i) = -(?!)(a;, i) < v{x, t)dx(l){x, t) + -dxv{x, t) \ (76) 

Multiplying ( I76] l by 2 (^ recovers the probability transport equation by the current 
velocity. 

5.3 Control equation 
Since 

dtf{x + dx, X, t) = dx dtv{x, t) + 0[dx^) {11) 

the continuum limit of ( I58ab calls for an expansion up to third order in dx: 

dx^ T (dtv) {x, t) = 

1 dx 

- [dxf [v (x, t) + — {d^v) {x, i)]2 2 dx d^ In 0(a:, i) 

1 da; 

+ - {dxf [v {x, t)- — (d^v) (x, t)f [1 - dx d. In (/.(x, t)] 

1 (jrn 

-- [dxf [v {x, t) + 3— (9^w) [x, t)Y [1 + dx d^ In (/)(j;, t)] + ©(dx)* (78) 

Unfolding the products we get into 

Tdtv{x,t)+v{x,t)d^v{x,t)^0 (79) 

which is the equation for the current velocity found in |,3j,4J . 

6 Examples: two and three-state systems 

To simplify the notation it is expedient to measure time in units of the traffic rate, 
tc = 1, and to consider a unit lattice spacing, dx = \. 
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6.1 Two-state system; comparison with optimal control for Glauber transition 
rates 

In the case of a two-state jump process, probability conservation yields a closed 
single equation for the evolution of the probability amplitude of the first state 

9*0(0, t) = - vT^WmF sinh ^^^ (80) 

We can, thus, proceed as in (TT\ and write the entropy production ( l60l) as a func- 
tional of the occupation probability 



of the first state. We obtain 



pit) = <p{o,ty (81) 



S^ dta{p,p) (82a) 



a{x, y) = ~y In ^^^i_^) ^ (82b) 

A straightforward and somewhat tedious calculation shows that the Euler-Lagrange 
equation specifying the trajectory of p corresponding to an extremal point of ( 
coincides with the one obtained by differentiating dSOb and using 



V{0,t) ,„^u VJO.t 



2 sinh ^™ tanh ^^ ( H^^ 



to write a closed expression for the evolution of the occupation probability of the 
first state 

v,p{l~p) +p'^ 

This is a consequence of the well-known fact (see e.g. }2l\) that any system of 
Euler-Lagrange equations is equivalent to a suitable deterministic control prob- 
lem via Hamilton-Jacobi theory. In writing ( 182b we fixed the reduced traffic to 
a constant, p-independent value as required by the optimal control equations of 
section |4] It is interesting to compare the entropy production in ( 182b with the one 
corresponding to a Glauber jump process with transition rate 

Ki(0|l) = ^^ (85) 

e" + 1 

since this is the jump process considered in [12|. We obtain 

/■*/ 
Sgi = I dtaGiip,p) (86a) 
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Fig. 1: Numerical illustration of the inequality {aGi—o-){x,y) > 0. We plot the dif- 
ference between the entropy production of the Glauber aci and of the Hamilton- 
Jacobi-Bellman optimal a jump processes for fixed value of the occupation prob- 
ability X = p as a function of y. The continuous curve corresponds to x = 0.05, 
the dotted to x = 0.35, the dashed to a; = 0.65 and the dash-dotted to x = 0.95. 
The inequality is proved analytically in the main text. 



'^Gi(x,y) = -yln— — — ■ — r > 



(86b) 



(1 -x){x + y) 

A preliminary observation is that whilst (I82t is well-defined for any < x < 
1, (|86] | restricts admissible controls to those satisfying the additional condition 
—X < y < I — X. Furthermore, the difference 



, N^ N , jx + y)[y^ + 2x{l - x) -y ^y4.x{l - x) + y^] ^^^^ 
[aci - <y)[x, y) ^y\n TT-TT. ' "1 (87) 



2a;2(l 



y) 



is positive-definite for any < a; < 1 and admissible value of y. To prove the 
claim it is sufficient show that 



9{x,y) 



{x + y)[y'^ + 2x{l- x) -y A/4a;(l - x) 



y" 



2a;2(i 



y) 



> 1 



(88) 



if y is positive while < g{x, y) < 1 if y is negative. Indeed ( 1881) reduces to 

y{{y^ + xi2 + y))-{x + y)^'ix{l-x)+y^}>0 (89) 

the argument of the curly brackets being positive whenever 

l-x-y+{x + yf >0 (90) 

The Glauber transition rate (185 1 corresponds to fixing the reduced traffic to 

1 



G(0,1) = 



2 cosh ■ 



(91) 



and the identification A(0, 1) — ~h. These modeling choices can be advocated 
with convincing physical arguments. The existence of "local" optimal controls for 
special choices of the reduced traffic suggests that the reduced traffic can be con- 
sistently, and perhaps should be, more appropriately thought as a configurational 
rather than a control parameter of a physical system. 
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6.2 Optimal control of the three-state system 

One of the simplest non-trivial application of the full-fledged optimal transport 
equations ( 1581) is to a system with three states. The probability amplitude is gov- 
emed by the equations 



dt(f>{0,t) = -(t){l,t) sinh 



1/(0, t) 



(92a) 



dt(j}{2,t) = <^(l,t) sinh 
while the current velocity satisfies 



V{l,t) 



dtViO,t) 



2 sinh 



vjo.t) 



tanh 



2 



2 — tanh 



2 V(0,t) 



a.t) 



2 sinh ^^ tanh ^^ 0(2, t) 



(92b) 



0(0, t) 



0(0, t) 0(1,0 



2 — tanh 



2 V(0.t) 
2 



0(1, t) 



= 



(93a) 



dtVil,t) 



2 sinh 



V(l:t) 



tanh 



vjht) 

2 



2 — tanh 



2 V(l,f) 



0(2, i) 0(1, i) 



2sinhMtanhM^(0,t) 



2 — tanh 



2 y(i,f) 

2 



0(1,0 



= 



with 



(1,0 = ^l-02(O,i)_02(2,i) 



0(1,0 0(2,0 J 



(93b) 



(94) 



The equations specify the minimal heat release protocol once we specify the initial 
and final states through ( [59] l. 

6.2.1 Numerical solution 



We considered the problem of minimizing the heat release during the transition 
between the states 



[m(0, 0), m(l, 0), m(2, 0)] = [0.9, 0.05, 0.05] ^ 
[to(0, 1), m(l, 1), m(2, 1)] = [0.05, 0.05, 0.9] 



(95) 



We integrated numerically the system governing the evolution of the occupa- 
tion probabilities [m(0, t), to(1, 0] and the discrete current [T/(0, t), T/(l,i)] us- 
ing Wolfram's Mathematica (editions 7.0 and 8.0.1) 1147 II default shooting method 
with at most 4 x 10^ number of iterations. We used an initial Ansatz for the pair 
V{Q, 0), V{1, 0) which we afterwards improved using the candidate solutions pro- 
duced by the shooting algorithm. The improvement of the Ansatz shifted further 
in time the blow up of the candidate solutions. After few manual iterations we 
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ria t ii> 



(a) Probability 



(b) Current velocity 



Fig. 2: Numerical solution of the three state process versus time. On the left (l2al) 
the behavior of the occupation probability of the first m(0, t) and the second 
m{l,t) state. On the right (l2b] i the current velocity components V{0,t), V{l,t) 
driving the transition. 



obtained the pair [K(0, 0), K(l, 0)] = [3.150, 1.156] which yields a smooth so- 
lution over the full control horizon [0, 1] matching the boundary conditions, see 
fig 121 The entire procedure takes few minutes using an old Pentium M powered 
laptop. We also checked that the same results can be recovered using the proba- 
bility amplitude equations but we noticed that this approach seems to increase the 
numerical stiffness of the problem. 

The probability of the first state is always decreasing while the probability of 
the intermediate state has a maximum at i = 0.5 and then symmetrically decays 
to match its boundary condition. The constant line corresponds to 0.05 and is 
plotted to emphasize the convergence of the two occupation probabilities to the 
the assigned boundary condition. The optimal controls are fast varying near the 
boundaries while varying more slowly in the bulk of the control horizon. This 
is the characteristic behavior of optimal controls of thermodynamics functionals 
described in Il41lll2l i3..4..2j. Since the Gibbs entropy vanishes for this symmetric 
transition the heat release coincides with the entropy production. The numerical 
evaluation of dSOl l for the three state system then yields 



/3Q* =S,= 4.096 



(96) 



7 Conclusions 



We have shown that if the traffic is constrained to a constant value the driving 
function governing the transition between two states with minimal heat release in 
a finite time horizon obeys Hamilton- Jacobi-Bellman type equations. The corre- 
sponding dynamics is short-ranged as it permits only jumps to nearest-neighbor 
sites. Our results are the counter-part for Markov jump processes of the entropy 
bounds derived in |2| for the Langevin dynamics. In this sense, they can also be 
regarded as refinement of second law of thermodynamics for Markov jump pro- 
cesses. We have also shown that entropy production is not convex in the traffic. 
As a consequence, interpreting traffic as a control parameter results into singular 
optimal control strategy. This means that any eventual heat release occurs in a 
zero Lebesgue measure time interval. It is worth noticing here that jump protocols 
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solutions survive also in the Langevin continuum limit. Namely, if we express as 
in JU the heat functional in terms of the current velocity it is readily seen that 
the lower bound provided by the variation of the Gibbs-Shannon entropy becomes 
tight for a jump protocol for which the current velocity is vanishing throughout 
the control horizon. This latter condition then just means that the drift must re- 
main in equilibrium so that the resulting protocol is precisely the one found here 
in the case of unconstrained traffic. In the Langevin case, this kind of solutions 
for the heat functional can be ruled out a-priori by restricting admissible control 
to smooth diffusions. The current velocity representation also evinces the subtle 
problems in which we incur when turning to thermodynamic work minimization. 
For the work it is not possible to find a minimizer by restricting admissible pro- 
tocols to smooth diffusions because of the non-coercive functional dependence 
on the control occasioned by the boundary cost term specified by the internal en- 
ergy variation. Justifying in which sense variations for fixed terminal value of the 
intemal energy may still define to an optimality criterion for a cost functional 
compatible with the interpretation of non-equilibrium thermodynamic work calls 
then for a more refined analysis such as that carried out in \^ . An comprehensive 
discussion of these results can be found in 



To summarize, we can say that the analysis of optimal control of Markov jump 
process helps to shed further light on the interpretation of the result previously 
available for the Langevin dynamics. Restricting admissible controls to those de- 
termining the driving function provide then a simplified model for the continuous 
limit Langevin dynamics. Such an approach is particularly beneficial because, as 
exemplified by the three state jump process, numerical solutions of the Hamilton- 
Jacobi-Bellman equations become extremely computing time inexpensive in com- 
parison to the Monge-Ampere-Kantorovich scheme solving the Langevin dynam- 
ics IJl- Still, the three state Markov jump process captures relevant qualitative 
features of the Langevin dynamics. The discrete dynamics is therefore very use- 
ful for capturing the qualitative features of the optimal control of non-equilibrium 
thermodynamics statistical indicators, such as higher moments of the heat func- 
tional itself, the Hamilton-Jacobi-Bellman equations thereof are not amenable, at 
least as far as it is currently known, to any fast integration scheme in the Langevin 
limit. 
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Appendices 

A Mean backward derivative 

Local detailed balance ^ ensures that the chain of equalities 

E„ , , I _ , - m(x, t) , 

^^^P.,,,4x,t + dt|x,t)x-^^^-l^ (97) 

Expanding in Taylor series, neglecting terms 0{dt^) and using the Master equations {4), (|5j 
yields the proof of l|8j. 

B Evaluation of tlie Kullbacli-Leibler divergence (TSi 

The first term on the right hand side of ( 115) is vanishing on average 



J to ,„c. 



xSS 



x,yS5 

The second addend 



Kt(x|y)-Kt(y|x) 



m(x, f) 



m(x, t) = (98) 



t^O '^^ (?*!?'-) ^*o ,% K,(y|x)m(x,i) 



carries two contributions. The first is the Shannon-Gibbs entropy of the transformation 
n 

x.ySS 

f ' dtY, lnm(y,i)[Kt(x|y)m(y,i)-Kt(ylx)m(x,f)] 



r'^^El-^K.(x|y)m(y,t) 



x.yeS 



m/(Ctj) 



= - / ^ dt y" In m(y, t) ftrnfy, t) = -E« In ^ ;!''\^ (100) 

Jto fts mo(CtJ 

by probability conservation eq. l |19) . The second specifies the thermodynamic heat and is given 
in the main text. 



C Alternative treatment of the variational problem 

We may state the optimization problem directly for the heat functional using as independent 
control the reduced traffic and the driving functional A: 

PQ^ dt' ^ A(x,5,t)G(x,i,t)e^^^m(5,t) (101) 

''*° x,SGS 
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The corresponding value function in such a case is the solution of the backward Kolmogorov 
equation 

[(at + L)J](X,t) + ^A(x,5,f)G(x,S,t)e^^^ =0 (102) 

xgS 

where, as in the main text, we require the variation of J to vanish at tf as we regard J(-, t/) as 
afunctional of m(-,i/). By Dynkin formula (see e.g. |27|) we must have 

I3Q = J2 {Mx, to) J{i, to) - m{i:, tj) J(5, tf)} (103) 

X 

whence using the boundary conditions ti follows that 

^Q' = ^m(i,t„) j'(X,to) (104) 

X 

The weak sense variation of i ll 02) yields the conditions 

= J(x,f) - J(5,i) + A(x,5,f) (105a) 

C A(x,x,t) A(x.ii,t) 1 J(x,f) - J(x,f) 

0=<m(x,t)e 2 +m(x,f)e ^ ( % — ^^ 

f A(,,x.t) A(x,x,t) 

+ <e 2 m(x, i)— e 2 m[x,t) 

A(x X t) ( A(x.x.t) A(x,x.t) ) 

+ ^ ' ' ^ <^ e^5~ m(5, f) + e ^5~ m(x, f) \ (105b) 

Upon inserting the Ansatz 

A(x, X, i) = F(x, X, t) + In m(x, i) — In 7Ti(x, t) (106a) 

J{x,t) = B{x,t) -lnm{x,t) (106b) 

we recover the stationary point equations ( 152) . 
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